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Numerical modeling of dependences "critical current - external magnetic field" for geometrically 
symmetric — tv Josephson junctions is performed. The calculation of critical current is reduced to 
non-linear eigenvalue problem. The critical curve of the contact is obtained as an envelope of the 
bifurcation curves of different distributions of the magnetic flux. 

The structure of vortices in contact is observed explicitly and the dependence of the basic physical 
characteristics of these vortices on junction's length is explored. 

The comparison of numerical results and known experimental data shows good qualitative and 
quantitative conformity. 
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I. INTRODUCTION 

The opportunity of existence and interest of studying 
of — 7r Josephson contacts (further on referred as 0-7rJ J 
s) have been shown in classical papers 0, 1 (see the 
short history of the problem in 3}). For last years in 
the specified area significant number of theoretical and 
experimental works are published (see for example, - 
[10| and references therein). 

In this article a numerical modeling of some experi- 
mental data [8| concerning the relations "critical current 
- external magnetic field" for such a junctions, is carried 
out. 

For this purpose the transitions from Josephson to re- 
sistive regime by change the external current 7 are math- 
ematically interpreted as a bifurcation of one of possible 
static distributions of the magnetic flux ip(x) existing in 
the junction for given external magnetic field h e (from 
now on all the variables are dimcnsionlcss). 

Every solution of non-linear boundary value problem 
(BVP) for <p(x) generates some regular Sturm-Liouville 
problem (SLP). The minimal eigenvalue (EV) Ao of SLP 
allows to decide whether the concrete solution is stable 
or unstable in linear approximation. Note, that such 
approach for the traditional model of "one-dimensional" 
long josephson junctions (JJs) was proposed in (Til. H3|. 

In such approach the critical current 7 C of some fixed 
distribution is that value of 7, for which the minimal 
EV Ao is vanishing at the given external magnetic field 
h e . In order to calculate points of bifurcation 7 C a non- 
linear eigenvalue problem is formulated, such that the 
current 7 is considered as an EV under the rest given 
junction's parameters including Ao [Hj|, Since for 

given external magnetic field the static BVP can have 
more than one solution, then the critical current of JJ 
is determined as maximal of critical currents of possible 



stable distributions. Thus the critical curve (CC) of the 
junction can be constructed as an envelope of bifurcation 
curves (BCs) of different distributions of the magnetic 
flux. 



II. STATEMENT OF THE PROBLEM 

Let consider 0—tt junction of length 21, 1 < 00, disposed 
along axis x, such that 0-junction (0JJ) is to the left and 
7r-junction (7rJJ) is to the right of the barrier 1 = ( e 
(—1,1). The junction is geometrically symmetric if £ = 0. 
In case of overlap geometry the non-linear BVP for static 
distributions of magnetic flux <p(x) can be written as 



-fxx + jc(x) simp - 7 = 0, 

xe H,c)u(c,0, 

V>(C-e) = <p(C + e), 
tp x {(-e) = tp x (( + e), 
ip x {±l) -h e = 0,, 



(la) 

(lb) 
(1c) 
(Id) 



where jc (x) is the amplitude of Josephson current 

jc(x) ■■ 
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Possible solutions of ([T]) belong to a class C : [— I, I], i.e., 
continuously differentiable (smooth) functions. Physi- 
cally this means that the magnetic field inside the junc- 
tion is continuous only, and at the sewing point x = C, 
continuity conditions (|lb[) and (fTc)) are fulfilled [H-Q- 

One can check that for 7 = the common solution 
of equations (|la[) in the corresponding intervals is ex- 
pressed by the elliptic functions, see [2(| ■ Under the given 
model's parameters p — {I, h e , 7} the integration con- 
stants can be obtained by means of two boundary con- 
ditions (jldp and two continuity conditions (|lb[) and (fTc|) 

For small magnitudes of the current I7I << 1 the so- 
lution of |T]) can be derived applying perturbation the- 
ory methods. But in all cases the solution of arising 



2 



non-linear algebraic systems can be obtained only nu- 
merically. Such approach for deriving bifurcation curves, 
which correspond to variation of parameters p, even 
7 = 0, meets serious difficulties. Therefore, in the present 
paper the direct numerical modeling of the problem ([1]) 
is preferred. 



III. STABILITY OF STATIC SOLUTIONS 

Further we shall suppose continuous dependence of so- 
lutions ip(x,p) of ([T]) on parameters p 6 V, V C K 4 . By 
now, the solutions dependence on p we shall mark only 
if it is necessary. 

In order to analyze the stability of solutions of (fT|) 
under parameters variation we put in correspondence to 
each solution ip(x,p) a regular Sturm-Liouville problem 



ditions for the total energy functional of the junction 



-ipxx + q{x)4> = Xtp, 

xe H,C)u(C,0, 
^(C-e)-V(C + e)=0, 

lMC-e)-Vfc(C + e) = 0, 

V»*(±0 = o, 

l 

i} 2 {x)dx -1 = 0, 

-l 



(2a) 

(2b) 
(2c) 
(2d) 

(2e) 



where the potential q(x) = jc{x) cos <p(x). Here A is a 
spectral parameter and (|2e() represents the norm condi- 
tion. 

It is easy to prove the existence of non-degenerate dis- 
crete bounded below spectrum {A„}^° of problem @, see 
[li ]. Then the bifurcation equation of order n of some 
solution ip(x,p) is 



A n (p) = 0, n = 0,l,. 



(3) 



It is convenient to interpret the last equality geometri- 
cally as a surface in the parameters' space V . Each point 
on this surface is a bifurcation (critical) point of order 
n for some solution of ((!]) . For fixed value of two of pa- 
rameters, the equation ([3]) describes a curve on the plane 
defined by other two parameters — a bifurcation curve 
(BC) of order n of considered solution. 

In this paper we mainly consider bifurcations of order 
of solutions of (JTJ) by change basic physical parameters — 
the external magnetic field h e and the external current 7. 
An implicit equation of BC for some solution on (h e ,j) 
plane is 



A (ft e , 7) = 0. 



(4) 



For OJJs bifurcations of static distributions of mag- 
netic flux under variation of the length / were studied, 
for example, in articles [Hj], [l7| . 

From mathematical point of view, the problem (fTJ) can 
be considered as a collection of necessary extremum con- 



F[<p] 



\<Px + ic{x) (1 - coacp) - 7^ 



dx — h e Aip 

on a set of smooth on [—1, 1} functions <f(x). 

Especially, the equations (|la|) are Euler-Lagrange 
equations for ([5]) at corresponding subintcrvals, and 
Weierstrass-Erdmann conditions yield (|ld[) , (|lb[) . and 
(jTc[) , see [l9j]. The total magnetic flux through the junc- 
tion Aip is defined traditionally 



Aip = / ip x dx = ip(l) — (p(— I). 



(6) 



-/ 



SLP © has to satisfy sufficient extremum conditions 
for ([5|), see [l9|. If Ao > for some solution ip(x) of 
([T]), then for the same solution the second variation of 
F[ip] is positive definite and the functional ([5]) possesses 
minimum. Therefore the solution ip(x) is stable in linear 
approximation. 

If Ao(p) = for some p e V then the solution <p(x,p) 
has at p a bifurcation of zero order, and the equations 
(|2a|) represents the Jacobi equations for the functional 
(J5j) . For brevity we shall call solutions at any parameter's 
bifurcation point as B-solutions. 

In this article we apply the method for numerical con- 
struction of relations (f?)) proposed in [IH; [l4| (see also 
review [Hi). At given A both equations ([1]) and ((2|) are 
considered as unique non-linear eigenvalue problem with 
spectral parameter h e or 7. To find eigenvalues and 
eigenf unctions ((p(x),ip(x)) we use an algorithm based 
on continuous analogue to Newton method [22]. The 
discretization of corresponding linearized BVPs is per- 
formed by means of spline-collocation method. Such 
approach leads to difference scheme with block three- 
diagonal matrix. The scheme's accuracy assessment, de- 
termined by Runge's method on the sequence of uniform 
meshes with steps h, h/2 and h/4, is 0(h 4 ). 



IV. DISCUSSION OF NUMERICAL RESULTS 

A. Non-bifurcation static solutions 

First let's consider some solutions of (p}, which exist 
far enough from the bifurcation point of the magnetic 
field h e or/and current 7, and behavior of their main 
physical characteristics when parameters tend to the crit- 
ical values. 

On Fig. [T] the main stable states (of minimal energy) 
of magnetic flux tp(x) for 0-7rJJ at two values of the junc- 
tion's length 21 = 2 and 21 = 7 in the field h e = and 
under the current 7 = are demonstrated. Further, for 
convenience, we shall use notation $ ±1 for mentioned 
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FIG. 1: Magnetic flux in 0-ttJJ 
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FIG. 2: Magnetic field in 0-ttJJ 



above solutions, where $ x is a fluxon and $ 1 - anti- 
fluxon. 

Corresponding distributions of magnetic field (deriva- 
tives tp x (x)) are shown on Fig. [21 As in the case of ring 
7r-contact, see [Ullj], the main states have finite total mag- 
netic field (A</? » ±0.147 at 21 = 2, and Aip w ±0.468 for 
21 = 7, respectively), finite energy (i 7 '[<£ ,±1 ] « —0.04, and 
F[<f> ±1 ] w —0.583) and contain a half-integer numbers 
jV[$ ±:l ] = 0.5 of fluxons (magnetic flux quanta), where 
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FIG. 3: Bound states in 0-7rJJ at h e = 2 



-3 -1.5 1.5 3 

Distance 



FIG. 4: Bound states in 0-7rJJ at h e = 3.8 
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FIG. 5: Dependence \o{h e ) for 0-7rJJ of length Z = 7 at 7 = 

m 

1 

N[<p] = ±J<p(x)dx, (7) 

-/ 

is the average magnetic flux trough the junction. 

At the center x = we have 4> ±1 (0) = ±7r/2 (see Fig. 
[IJ. Moreover, the conditions of geometrical symmetry 
are fulfilled ^(x) = -^(x), $^ 1 (x) = -®l(x), and 

Except vortices demonstrated on Fig. (TJ at h e = 
and 7 = there exist Meissner's states Mo (a;) = and 
M v (x) = 7r (+2nk, k — 0, ±1, ...). Such distributions 
have zero total magnetic flux Aip and zero energy 
but N[M ] = and N[M^] = 1. Note, that the po- 
tential of SLP induced by M , is q{x) — jc( x ) an d 
q(x) — —jc(x) for M„. Hence no one of Meissner's so- 
lutions can be stable in contrast to traditional 0JJ and 
7rJJ, where q(x) = 1 for Mq (stable) and q(x) = — 1 for 
Mjr (unstable). In particular Ao « —0.3 at 21 = 2 and 
Ao ~ —0.86 for 21 — 7. Numerical experiment confirms, 
that Meissner's solutions remain unstable in all admis- 
sible range of magnitudes h e and 7. Thereby in case of 
0-7rJJ the traditional Meissner's branches are missed on 
CCs, as it is observed in the experiment §]. 

Increasing value of h e arise new more complicated 
bound states of the magnetic flux in the junction — pure 
and mixed chains of fluxons and antifluxons. Pure chains 
of vortices are "composed" only of fluxons or antiflux- 
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FIG. 6: Dependence Ao(7) for 0-7rJJ of length I — 7 at h e — 
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FIG. 7: Magnetic flux distributions ip{x) at 7 = 



ons. Mixed chains are obtained after non-linear inter- 
action between fluxons and antifluxons. To mixed so- 
lutions of type $ m $ _m correspond symmetric <J> _m <I> m , 
to = 1,2,.... On Fig. [3] chains of vortices of internal 
magnetic field ip x (x) for the field h e — and at the 
current 7 = are demonstrated: one-soliton distribu- 
tion $ , system of three non-linear ly interacting solitons 
(7V[$ 3 ] = 2.5), and also $ 1 <i>~ 1 , pair soliton - anti- 
soliton, with TV^ 1 ^ 1 ] « 0.005; the symmetrical pair 
yields JVf* -1 * 1 ] » 0.995. Similarly, on Fig. H 
a pure chain of fifth solitons <E> 5 (7V[<I> 5 ] = 4.5), and 
mixed chain $ 2 <i> -2 of two solitons and two antisolitons 
(N[$ 2 $~ 2 } fn 4.002 at h e = 3.6), pinning at the inhomo- 
geneity at the point ( = are shown. For the symmetri- 
cal pair $- 2 $ 2 we have 7V[$~ 2 $ 2 ] w 2.998. 

On Fig. [5] the dependence of minimal EV of ([2|) on 
the external magnetic field h e for the junction of length 
21 = 7 at 7 = is represented. In accordance with (|4|) 
the zeroes of curves Xo(h ei 0) are bifurcation points of 
the corresponding solutions. Each curve possesses two 
zeros which correspond to the lower h m i n and upper h max 
critical magnetic field for concrete distribution ip(x). The 
distance Ah = h max — h m i n represents the domain of 
existence of the distribution ip(x) on the field h e at 7 = 
0. In particular, for ^-distribution calculations yield 
hmin ~ -1.8 and h max « 2.2, thus Ah « 4. 

Similarly, on Fig. [5] the behavior of the minimal EV 
of SLP by change the current 7 for 21 — 7 and h e = 
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FIG. 8: Magnetic field distributions ip x (x) at 7 = 
" 0-7IJJ, <*>': 21 = 7, he = 
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FIG. 9: Magnetic flux distributions tp(x) at h e — 



is demonstrated. Each upper-down curve corresponds to 
two symmetrical vortex distributions with the same Ao 
and same energy F, but with opposite signs of number 
of fluxons and Atp. The curve in the upper half-plane 
corresponds to single fluxon distributions and 
Points Bo (I7I k 0.6) are points of bifurcation of men- 
tioned solutions varying parameter 7. At these points 
arises confluence of $ 1 with Meissner's distribution Mo, 
which is much deformed by the current. For the given 
junction, the distance A7 w 1.21 between points Bq de- 
fines domain of stability of distributions $ x via current 
7 in the field h e — 0. 

At points B\ the first EV Ai of SLP is vanishing. This 
means that there is a bifurcation of unstable distribution 
to another unstable one. In the experiments it is difficult 
to discover such transitions because of small life time [T3] • 

Fig.[7]illustrated the deformation of the basic fluxon $ : 
by external field in 0-7rJJ of length 21 = 7 at the current 
7 = 0. Corresponding distributions of magnetic field 
(p x (x) in contact are given at Fig. [51 It is well observed, 
that the main changes are localized at the neighborhood 
of contact's ends. The center ip(0) of fluxon is conserved. 

Similarly, on Fig.[5]the deformation of fluxon by the 
current 7 at zero field h e is represented. Positive current 
7 displaces the fluxon's graph upwards (the maximum of 
corresponding soliton <p x (x) moves to the left), and the 
negative vice versa — the graph of f(x) moves down, and 
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FIG. 10: BCs of the main fluxon $ 1 

the graph of soliton <p x (x) — to the right. 

The basic numerical characteristics of the first several 
non-bifurcation magnetic flux vortices in contact at 21 = 
7 and 7 = are given in the Table |TJ 



Type 


h e 


Ao 


N[<p] 


A<^/2tt 


V(0)/Jr 


FM/8 


$-1 





0.772 


-0.5 


-0.468 


-0.5 


-0.586 


$i 





0.772 


0.5 


0.468 


0.5 


-0.586 




1.61 


0.0275 


0.885 


1.7 


0.664 


-1.349 




1.61 


0.0267 


2.115 


1.7 


2.336 


-1.346 


(J) 3 


2 


0.108 


2.5 


2.582 


2.5 


-1.883 




3.41 


0.015 


2.897 


3.773 


2.771 


-5.144 




3.41 


0.012 


4.103 


3.773 


4.229 


-5.140 


$ 5 


4.6 


0.111 


4.5 


5.049 


4.5 


-9.378 




5.28 


0.029 


4.963 


5.882 


4.866 


-12.22 


$ 3 $~ 3 


5.28 


0.023 


6.037 


5.882 


6.134 


-12.217 


$ 7 


6.4 


0.081 


6.5 


7.072 


6.5 


-18.00 



TABLE I: Magnetic flux vortices in 0-7rJJ 



One can see that the pure magnetic flux distributions 
match to half-integer values of number of fluxons ([7]) 

JV[$ n ]=nT-i 

where the minus sign correspond to n > and plus — to 
n < 0, while for mixed fluxon states we have 

7V[$«d>-«] + AT[$-"$ n ] = 2n T- -■ 

Let's note that the quotients <p(0)/27r fulfil the similar 
relationships as well. 



B. Bifurcations of static solutions 

In this article we shall mainly study bifurcations of 
static solutions varying external field h e and external cur- 
rent 7. 

The locus on the plane (/i e ,7) satisfying equation ([J), 
we shall call a bifurcation curve. A unique BC corre- 
sponds to each solution of non-linear BVP ((T|). 
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FIG. 11: BCs of pairs 
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FIG. 12: BCs of distributions in 0-7rJJ contact of length 21 = 7 



For example, on Fig. [TU] the BC of the vortex <I> 1 in 
0-7TJJ of length 21 = 7 at 7 > (solid line) and 7 < 
(dashed line) is demonstrated. The distance Ah e w 4 be- 
tween zeroes h m i n w —1.8 and h max k, 2.2 represents the 
region of existence of $ x by change the external magnetic 
field h e . The graphs of derivatives (p x (x) of solutions in a 
small neighborhood of bifurcation points h m i n and h max 
are shown on Fig.[5]by dotted and dash-dotted lines cor- 
respondingly. 

Similarly, on Fig. QT] the BCs for the pair C^ - ^ 1 , 
are demonstrated. The existence domain of dis- 
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FIG. 13: CC for 0-ttJJ of length I = 7 
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tributions via field Ah = h max — h. m i n is separated by 
intermediate points hi and h 2 , for which 7 cr = 0, on 
three subdomains. 




Symmetric O—Ji contact 
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FIG. 14: Disposition of the point An(0, 7) as a function of 
length I 

Let us note, that for h e e (h m i n ,hi) and h e £ 
(^2, h max ) the corresponding domains of existence of dis- 
tributions via current 7 are lied over/under the horizon- 
tal axes 7 cr = 0. Thus when |7| increases from zero, the 
mixed vortices at the mentioned domains have a current 
of "birth" (for given h e these are points on dotted curves) 
and a current of "destruction" — points on solid curves. 
For h e = h m i n and h e — h max the currents of "birth" 
and "destruction" coincide. 

At the points h e = hi and h e = h 2 the critical current 
of distributions and (I) 1 ^ -1 is equal to zero. 

The effect under consideration presents also in JJ with 
micro- inhomogeneities of resistive or shunt type [l5j |. 
Maybe, this effect can be determined experimentally us- 
ing methods proposed in [2lj . 




External magnetic field 



Because of the symmetry, only the right upper quarter of 
CC is demonstrated. Let us note, that CC in symmetric 
0-7rJJ represents a set of branches, which sequentially 
correspond to pure and symmetric mixed distributions. 

An important peculiarity is the shift of the critical cur- 
rent maximum 7 cr to the right for h e > (to the left 
for h e < 0) from the vertical line h e = 0. The point 
70 = "far (he = 0) essentially depends on the contact's pa- 
rameters — (half)length I and location £ of the barrier. 

For the geometrically symmetric junction (£ = 0) the 
dependence 70 (0 is given on Fig. Q31 For small I the 
current 70 non-linearly quickly enough tends to zero. At 
large length I the current 70 asymptotically diverges to 
the constant value 700 ~ 0.637, while the total magnetic 

flux Atf(2l) -> 7T. 




0-JtJJ: 21 = 7 
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FIG. 16: Bias of intersection points of $ & $~ 1 critical 
curves for different C, 

A comparison of results of numerical modeling of CC 
(entire curve) with available experimental data [1[ for 
short 0-7rJJ of length 21 — 1.3 is demonstrated on Fig. 
IT5l Note, that for such length the contribution of BC of 
mixed states on junction's CC can be neglected. 

A good qualitative correspondence between numerical 
and existing physical experiments is demonstrated. A 
plausible reason for some quantitative differences can be 
the existence of different in-symmetries in experimental 
samples, as £ ^ and etc. For example, at Fig. [16] the 
influence of the shift ( ^ on BCs in case 21 = 7 
is demonstrated. For £ = the intersection point Cq 
correspond to /i e |c=o = 0; while the inhomogeneity bias 
to the left or right (£ = =f0.5) leads to bias of corre- 
sponding intersection points to the right or left 
(ft e |c=T0.5 ~ ±0.17) correspondingly. 



FIG. 15: CC for 0-vrJJ of length I = 1.3 

All possible BCs of distributions of magnetic flux in 
0-7rJJ of length 21 — 7, existing at widely differing h e , 
are given on Fig. The critical curve of the junction 
contains all points of BCs with maximal modulo critical 
current at a given field h e , i.e., the contact's CC is derived 
as envelope of BCs of all distributions. On Fig.[T3]CC of 
the contact mentioned above is illustrated for h e < 10. 



V. CONCLUSIONS 

We have numerically modelled experimental data Q 
for symmetric — it Josephson junctions. 

Essentially, some stable bound states of magnetic 
fluxes and corresponding magnetic fields distributions are 
derived and demonstrated. For 0-7r JJ of length 21 — 7 we 
obtain and illustrate relations of minimal eigenvalue Aq 
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on external magnetic field h e and external current 7. We 
illustrate numerically the important role of average mag- 
netic field in the JJ as a measure of fluxon's number in 
a vortex chains. Special BCs and CCs for some fluxons 
and antifluxons are represented, too. 

For relatively short 0-7rJJ a good qualitative confor- 
mity between numerical and existing physical results 
takes place. We propose that some quantitative differ- 
ences can be caused on different in-symmetries in exper- 



imental samples. 
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